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ABSTRACT 

Methods  for  studying  the  bifurcation  behavior  of  tubular  reactors  have 
been  developed.  This  involves  the  application  of  static  and  Hopf  bifurcation 
theory  for  PDE's  and  the  very  precise  determination  of  steady  state 
profiles.  Practical  computational  methods  for  carrying  out  this  analysis  are 
discussed  in  some  detail.  For  the  special  case  of  a  first  order,  irreversible 
reaction  in  a  tubular  reactor  with  axial  dispersion,  the  bifurcation  behavior 
is  classified  and  summarized  in  parameter  space  plots.  In  particular  the 
influence  cf  the  Lewis  and  Peclet  numbers  is  investigated.  It  is  shown  that 
oscillations  due  to  interaction  of  dispersion  and  reaction  effects  should  not 
exist  in  fixed  bed  reactors  and  moreover,  should  only  occur  in  very  short 
"empty"  tubular  reactors.  The  parameter  study  not  only  brings  together 
previously  published  examples  of  multiple  and  periodic  solutions  but  also 
reveals  a  hitherto  undiscovered  wealth  of  bifurcation  structures.  Sixteen  of 
these  structures,  which  come  about  by  combinations  of  as  many  as  four 
bifurcations  to  multiple  steady  states  and  four  bifurcations  to  periodic 
solutions,  are  illustrated  with  numerical  examples.  Although  the  analysis  is 
based  on  the  pseudohomogeneous  axial  dispersion  model,  it  can  readily  be 
applied  to  other  reaction  diffusion  equations  such  as  the  general  two  phase 
models  for  fixed  bed  reactors. 
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SIGNIFICANCE  AND  EXPLANATION 


Although  the  theory  of  static  and  Hopf  bifurcation  for  nonlinear 
distributed  parameter  systems  has  been  essentially  developed,  computational 
procedures  easily  used  by  applied  mathematicians  and  engineers  are  required  to 
bring  these  results  into  practice.  The  object  of  the  present  study  is  to  show 
how  the  static  and  Hopf  bifurcation  behaviour  for  highly  nonlinear  processes 
may  be  determined  as  a  function  of  process  parameters.  Special  numerical 
methods  are  required  to  handle  the  very  stiff  steady  state  equations  often 
encountered.  Then  special  methods  are  demonstrated  for  mapping  the  parameter 
dependence  of  the  bifurcation  points.  The  direction  and  stability  nature  of 
the  oscillatory  solutions  arising  at  Hopf  bifurcation  points  is  not  determined 
here.  However  Poore  and  Heinemann  {1980]  describe  an  algorithm  for  this 
calculation. 

The  computational  procedures  are  illustrated  by  application  to  an 
important  problem  found  in  chemical  engineerings  the  dynamic  behaviour  of 
tubular  reactors.  Regions  in  parameter  space  are  determined  showing  where 
static  and  Hopf  bifurcations  (both  single  and  multiple)  occur.  The 
engineering  significance  of  the  results  for  this  example  problem  is  discussed. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


THE  BIFURCATION  BEHAVIOR  OF  TUBULAR  REACTORS 
Klavs  F.  Jensen^  and  W.  Harmon  Ray 

X .  INTRODUCTION 

In  fixed  bed  reactors  the  reactants  flow  in  gaseous  or 
liquid  form  through  a.  vessel  (generally  cylindrical  in  shape) 
packed  with  solid  catalyst  particles  within  which  the  reaction 
takes  place.  In  the  absence  of  catalyst  packing,  the  vessel 
only  serves  to  confine  the  reaction  medium  and  it  is  then 
known  as  the  "empty"  tubular  reactor.  Naturally,  the  complex 
transport  and  reaction  processes  allow  for  different  levels  of 
sophistication  in  the  mathematical  modelling  necessary  to  design 
and  optimize  fixed  bed  reactors  [1-5].  Since  these  reactors 
moreover  are  used  in  a  vast  number  of  industrially  important 
reactions  [1,  Table  1],  they  have  received  considerable  atten¬ 
tion  during  the  last  two  decades.  In  particular,  the  occurrence 
of  multiple  steady  s t ate s,  t. a ve 1 1 ing  waves,  and  oscillatory 
states  in  these  reactors  has  been  a  focal  point  in  many 
theoretical  and  experimental  investigations  (cf.  [5-12]  for  an 
overview).  These  phenomena  were  predicted  by  theoretical 
analyses  before  being  consciously  noted  in  experiments. 

Table  1  gives  a  listing  of  experimental  studies  where 
multiple  steady  states  ( M ) ,  so-called  wandering  profiles  (W) , 
and  oscillating  behavior  (0)  were  reported.  Following 
Schmitz  [9]  we  have  omitted  studies  where  the  reactor  was 
of  the  recirculating  type  which  behaves  more  like  a 
continuously  stirred  tank  reactor  (CSTR)  than  a  fixed  bed 
reactor.  The  publication  dates  of  the  entries  in  the.  table 
clearly  demonstrate  the  recent  growth  in  experimental  evidence 
for  bifurcation  phenomena  in  fixed  bed  reactors. 
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Re  f e  ronce 

Experimental  System 

Rema  rks 

14  . 

Hegedus  et 
1977  [26] 

al  ■  , 

Oxidation  of  CO  on  Pt/AljOj 
in  an  isothermal  fixed  bed 
reactor 

M1 

IS. 

Schleppy  and  Shah, 
1977  [27] 

NO  reduction  with  CO  over 
fiberglass  supported  Ru  in  a 
nonadiabatic  fixed  bed  reactor 

M 

16. 

Butakov  and 
Shkadi nsk i i 
1978  [28] 

1 

Liquid  phase  decomposition  of 
dini troxydi e t hy Ini t ramine  in 
acetic  anhydride  in  a  tubular 
re  actor 

0 

17. 

Hlavacex.  and 
Votruba,  1978  [12] 

Oxidation  of  CO  on  Cu0/A1203, 
Pd/Al^Oj  and  Pt/Al203  in  an 
adiabatic  reactor 

M 

18. 

Oh  et  al . , 
[29] 

1978 

Oxidation  of  CO  on  Pt/Al203 
in  an  isothermal  fixed  bed 
reactor 

Ml 

19. 

.Ukus  et  al 
1979  ['30] 

•  t 

Oxidation  of  CO  on  Pt/Al203  in 
an  adiabatic  fixed  bed  reactor 

M 

20. 

Oil  et  al .  , 
[31] 

1979 

Oxidation  of  CO  on  Pt/Al2o3  in 
an  isothermal  fixed  bed  reactor 

M1 

21. 

Sharma  and 
1979b  [32] 

Hughe  s  , 

Oxidation  of  CO  on  a  CuO-catalyst 
in  an  adiabatic  fixed  bed  reactor 

M 

22. 

Hlavacek  et 
■ °G0  [33]' 

al  .  , 

Oxidation  of  CO  on  Pt/Al203  in  a 
deactivated  fixed  bed  reactor 

M 

23. 

Kalthoff  and 
Vortmeyer,  1980 
[34] 

Oxidation  of  ethane  on  Pd/Al203 
in  a  nonadiabatic  fixed  bed 
reactor 

M 

24  . 

Paspek  and 
1980  [35] 

Varma  , 

Oxidation  of  ethylene  on  Pt/Al2Oj 
in  a  nonadiabatic  fixed  bed 
reactor 

M 

25. 

Puszynski  and 
Hlavacek, 1980 
[36] 

Oxidation  of  CO  on  Pt/Al203  in 
a  nonadiabatic  fixed  bed  reactor 

M  ,  W 

M  :  multiple  steady  states 
0  :  self-sustained  oscillations 

W  t  special  investigation  of  " wande r ing " -prof i 1 e  . 

1)  short  catalyst  bed  -  approximately  8  catalyst  layers 
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Because  of  the  complex  transport  and  reaction  processes 
in  fixed  bed  reactors,  multiple  steady  states  are  generated 
by  various  kinetic  and  physicochemical  mechanisms.  The 
observed  multiplicity  behavior  in  the  liquid  phase  reactions 
(entries  4  and  6)  has  been  shown  to  follow  the  predictions 
of  a  plug  flow  reactor  in  a  recycle  loop.  Naturally,  the 
multiplicity  ot  states  in  the  autothermal  reactors  (entries 
10  and  13)  can  be  attributed  to  the  feedback  of  heat  through 
the  preheating  loop.  In  the  studies  by  Hegedus,  Oh  and  their 
coworkers  (entries  14,  18,  and  20)  the  multiple  profiles  stem 
from  multiple  steady  states  of  the  individual  catalyst  particles. 
Because  the  bed  contains  relatively  few  particles,  the 
investigators  are  able  to  realize  a  number  of  stable  profiles. 
Paspek  and  Varma  (entry  24)  also  attribute  the  multiplicity 
behavior  to  the  individual  catalyst  particles  and  explain  the 
phenomenon  in  terms  of  interactions  between  the  reaction  and 
intraphase  transport  processes. 

In  the  cases  corresponding  to  the  remaining  entries  in 

Table  1,  dispersion  effects  seem  to  be  part  of  the  underlying 

mechanism.  Hlavacek,  Votruba  and  their  respective  coworkers 

(entries  8,  11,  12,  17,  and  22)  studied  extensively  the  effects 

of  reaction  conditions  on  the  multiplicity  behavior  in  C0- 

oxidation  and  found  that  the  phenomenon  disappeared  beyond 

a  critical  length  (corresponding  to  Pe  ~  180)  where  the 

m 

dispersion  effects  became  insignificant.  However,  they  also 
observed  three  stable  profiles  in  their  adiabatic  reactor 
contrary  to  the  theoretical  prediction  from  the  pseudo- 
homogeneous  dispersion  model  of  a  maximum  of  two  stable 
profiles.  On  the  other  hand,  Sharma  and  Hughes  (entry  21) 
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studying  the  same  catalytic  system  could  not  realize  more  than 
the  predicted  two  stable  profiles.  These  authors  found  that 
a  two  phase  axial  dispersion  model  was  required  to  accurately 
model  their  experimental  data,  whereas  Schleppy  and  Shah 
(entry  15)  showed  for  a  different  reaction  system  that  the 
pseudohomogeneous  axial  dispersion  model  sufficed  to  fit 
the  observed  ignition  and  quench  behavior.  Kalthoff  and 
Vortmeyer  (entry  23)  also  used  a  pseudohomogeneous  model  but 
found  it  necessary  to  include  the  radial  porosity  and  velocity 
distributions  in  order  to  quantitatively  model  the  observed 
multiplicity  behavior. 

Wicke  and  his  coworkers  (entries  2  and  3)  demonstrated  the 
existence  of  travelling  wave  fronts,  so-called  wandering  profiles, 
which  moved  slowly  (linear  velocity  -10  ^cm/sec)  and  with  small 
changes  in  shape  through  the  reactor  for  changes  in  the  gas 
velocity.  For  high  gas  velocities,  the  front  was  blown  out, 
while  for  low  velocities  the  front  moved  upstream  to  the  reactor 
inlet.  An  intermediate  velocity  stabilized  the  front  in  the 
middle  of  the  reactor.  The  travelling  fronts  are  characterized 
by  steep  concentration  and  temperature  gradients  and  most 
likely  represent  transitions  between  multiple  steady  state 
profiles.  Transients  reported  by  Sharma  and  Hughes  (entry  21) 
as  well  as  by  Puszynski  and  Hlavacek  (entry  25)  nicely 
demonstrate  the  moving  front  structures  which  come  about 
during  the  transition  between  states. 

Oscillatory  states  in  tubular  reactors  are  reported  in 
three  experiments,  entries  1,  9  and  16,  but  only  in  the  last 
case  are  the  oscillations  linked  to  a  mathematical  model, 
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namely  the  axial  dispersion  model.  The  lack  of  experimental 
evidence  of  oscillatory  profiles  in  fixed  bed  reactors  may  be 
explained  by  the  fact  that  the  characteristic  time  for  thermal 
transport  within  the  bed  is  so  much  larger  than  the  one  for 
material  transport,  i.e.  the  Lewis  number  is  so  large  that 
oscillations  are  not  expected  even  in  shallow  beds.  This  point 
will  be  treated  in  more  detail  below. 

The  question  of  uniqueness  of  the  solution  to  the  fixed 
bed  reactor  equations  has  been  considered  by  numerous  investi¬ 
gators,  notably  Amundson,  Hlavacek  and  their  respective 
coworkers.  The  axial  dispersion  model  with  a  single  first 
order  reaction  has  been  a  favorite  target  for  mathematical 
analyses  and  sufficient  conditions  for  uniqueness  have  been 
developed  by  applying  fixed  point  methods,  comparison  theorems, 
bifurcation  theory,  and  Liapunov  functionals.  These  contribu¬ 
tions  and  others  dealing  with  alternative  models  are  reviewed 
by  Jensen  [37],  Ray  [7],  Schmitz  19),  and  Varma  and  Aris  [11]. 
The  sufficient  conditions  for  uniqueness,  in  the  general  case 
of  a  nonadiabatic  reactor  with  unequal  Peclet  numbers  for  heat 
and  mass  dispersion,  show  (as  intuitively  expected)  that  the 
solution  will  be  unique  for  sufficiently  high  values  of  the 
Peclet  numbers,  large  heat  transfer  coefficients,  or  small 
values  of  the  Damkohler  number.  Extensive  calculations  by 
Hlavacek  and  his  coworkers  [38-43]  confirm  this  and  show,  in 
addition,  that  increasing  the  adiabatic  temperature  rise  or 
the  activation  energy  enlarges  the  region  of  multiplicity  and 
shilts  it  towards  lower  Damkohler  numbers,  while  an  increase 
in  the  reaction  order  reduces  the  region  of  multiplicity. 
Multiplicity  higher  than  three  is  possible  only  in  the  non- 


adiabatic  reactor,  where  five  steady  states  have  been 
calculated  140,44,45).  Recently,  Kapila  e_t  a_l_.  (46)  have 

shown  by  using  activation  energy  asymptotics  that  as  many  as 
seven  steady  states  may  exist  in  the  limit  of  large  activa¬ 
tion  energy. 

Also  the  travelling  reaction  fronts  have  generated  much 
theoretical  interest.  Vortmeyer  £t^  a _1.  [47-49]  consider  the 

reactor  infinitely  long  which  seems  to  be  a  reasonable  assump¬ 
tion  because  the  concentration  and  temperature  changes  in  the 
front  occur  over  very  small  distances  compared  to  the  reactor 
length.  Gilles  [50],  in  addition,  approximated  the  reaction 
rate  over  the  front  zone  by  a  Gaussian  distribution.  Rhee  e_t  a  1  . 
[51,52]  employed  two-phase  cell  and  continuum  models  and 
developed  explicit  formulae  for  the  velocity  similar  to  that 
for  a  shock  layer  in  a  nonreactive  system.  The  above 
approaches  all  showed  good  agreement  with  experimental  data. 

The  stability  of  the  steady  state  has  been  studied 
numerically  by  linearization  where  the  dominant  eigenvalues 
were  determined  from  either  a  collocation  or  Galerkin  approxi¬ 
mation  to  the  linearized  equations  [44,45,53,54].  Alternatively 
sufficient  conditions  for  stability  have  been  derived  through 
the  use  of  comparison  theorems  and  Liapunov  functionals  [55-61]. 

As  expected  from  the  physical  situation  and  formally  shown 
by  singular  perturbation  theory  [64,65],  the  pseudohomogeneous 
axial  dispersion  model  reduces  to  the  CSTR  model  as  the  Peclet 
number  becomes  very  small.  Therefore,  based  on  the  dynamic 
behavior  of  the  CSTR  [cf.  66,67]  one  expects  oscillations  to 
exist  in  short  reactors.  The  existence  of  such  oscillations 


in  the  ''empty”  tubular  reactor  has  been  determined  computationally 


by  Hlavacek  and  Hofmann  1391  and  Varna  and  Amundson  145b)  . 
However,  examples  have  not  been  calculated  for  fixed  bed 
reactors.  A  detailed  parametric  study  has  not  yet  been 


performed  for  either  type  of  reactor;  thus,  in  this  paper 
we  shall  show  how  such  a  study  can  be  made  by  using  bifurca¬ 
tion  theory.  We  shall  be  especially  concerned  with  the  effect 
of  the  Lewis  number,  i.e.  the  ratio  of  physical  transport 
thermal  time  constant  to  physical  transport  material  time 
constant.  This  parameter,  which  is  unity  for  the  empty 
tubular  reactor  and  much  greater  than  unity  for  fixed  bed 
reactors,  has  been  shown  to  have  a  striking  influence  on 
the  dynamics  of  chemically  reacting  systems  (cf.  [67]  and 
re f erences wi thi n ) .  Therefore,  the  analysis  will  have  practical 
interest  in  revealing  if  limit  cycles  are  at  all  possible  even 
in  snort  fixed  bed  reactors.  Although  the  bifurcation  analysis, 
i.e.  the  study  of  multiple  steady  states  and  oscillatory 
behavior,  will  be  based  on  the  p se udohomoge ne ous  axial 
dispersion  model  with  first  order  Arrhenius  kinetics,  the 
general  approach  readily  applies  to  other  reaction-diffusion 
equations  with  complex  rate  expressions.  The  dispersion  model 
is  a  particularly  good  example  for  illustrating  the  techniques 
since  the  modelling  equations  form  a  relatively  simple  set  of 
nonlinear  parabolic  differential  equations  which  are  capable  of 
predicting  multiple  steady  states  and  limit  cycles.  Moreover, 
stiffness  problems  are  often  encountered  in  the  numerical 
solution  to  the  steady  state  equations  which  means  that  one 
has  to  devise  versatile  and  efficient  algorithms.  Finally,  the 
previously  published  examples  of  multiple  steady  states  te.g. 
38-41,45)  make  it  possible  to  check  the  algorithms. 
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2.  Mom:?,Li::.  ^.cations 

In  t  h e  oa bc  of  an  irreversible  first  order  reaction,  the 
equations  of  the  pse udohomogeneous  axial  dispersion  model  are 
two  coupled  nonlinear  parabolic  partial  differential  equations 
for  the  reactant  concentration  and  the  temperature  [2,11]. 

The  equations  are: 


ac(?.'  ,  t 1 )  _  a  c  ( 7.  •  ,  t  * )  _  v  3c  (z '  .  t 1 ) 

EP  3t '  ‘  L  .  *2  z 


dz 


d  z 


kQ<l-Cp)c(z'  ,t 1  )  exp  [  - F,^/RT  (z  '  ,f)  ] 


(1) 


3T (z  •  ,  t  * ) 

lEppfcPf  +  a-£P!pscPs!~rd — 


3T2  (z  ’  ,t  ■  ) 


3T (z  *  , t  '  ) 


-  P,C  ,v 


3z 


£  p  f  z  3  z  ' 


♦  (1-E  )  (-£■-.)  k„c(z'  ,t')exp[-Es/RT(z'  ,f)  ’  ,  t ' > -T  <z ’  , t ' )  1 

p  0  A  u  j.  w 

(2) 


with  initial  conditions 


c(z 1  ,0)  =  c  (z • )  ,  T(z  •  ,0)  =  T. 


i  n 


(z  •  ) 


(3) 


Here  we  use  Danckwerts'  boundary  conditions,  even  in  the  analysis 
of  the  transients. 


-D.  --4— :-t—  =  v  [c  (f)  -  c(0,t')J 
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(da) 
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(5a) 


3T  (  z  '  ,  t  '  ) 
8z  * 
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-y- 


(5b) 


In  the  above  equations  c  and  T  represent  the  concentration 


and  temperature, 
coefficients,  while  v 


D  k  and  k 


are  the  longitudinal  dispersion 
is  the  gas  velocity.  Subscripts  f,  s. 


and 


denote  the  fluid  phase,  solid  phase,  and  reactor  wall 


respectively.  and  d^_  are  the  bed  porosity  and  diameter, 

kg  and  represent  the  usual  Arrhenius  parameters.  is 

the  overall  heat  transfer  coefficient  between  the  reactor  and 
cooling  medium.  The  effect  of  radial  heat  dispersion  may  be 
included  in  by  making  a  one-point  collocation  approximation 

to  the  radial  temperature  profile  [80]. 

The  equations  are  made  dimensionless  by  defining: 


c'co 

1 

[t  -t  J 

W  0 

A  t  .  _A_ 

X 

M 

II 

O 

O 

•  *2  = 

l  To  1 

RTq 

t  X » 

2w 

l  To 

RT0  '  RTq 

t  '  V 

2 

z  ’ 

Pp  — 

v  8. 
z 

‘  €  * 

P 

,  2  - 

t 

t 

pei 

dl~ 

V 

Pe  -  ~ 

tpfCpf 

Da 

in 

-c  )k  e~Y 

P  0 

t-AH) c0Y 

2 

kL 

V 

z 

• 

B  - 

pfcPfTo  ' 

(6) 


4VJ 


8  = 


v  d  p  C 
z  r  f  pf 


Le  - 


E„PfC  - 
P  f  Pf 


+  (1- C  ) p  C 


P ,C  ,e 
f  pf  p 


s  ps 


We  then  obtain  the  following  set  of  equations: 
3x 

C  jt-  Lx  +  £(*) 


where  C  denotes  the  capacity  matrix 


1 


0 

Le 


(7) 


(8) 
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and  L  denotes  the  linear  differential  operator: 


-1 

*■ 

Pe1x  0 

ji 

i 

0 
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0 

0 

„  2 

0  Pe"1 

dz 
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— 

The  nonlinearity  f(x)  is: 


f(x) 


Da(l-x^) exp 


X+x2/y 


B  Da(l-x^)exp 


1+x2/y 


+  gx 


2  v 


In  this  formulation  the  boundary  conditions  take  the  form: 


gz  1 


8^(0, t) 


-  Pe, 
gz  2 


x(  ,t)  -  o 


B2x(l,t) 


II  0 


3z 


x(l,t)  »  o 


In  the  special  adiabatic  case  with  Le  =  1,  the  modelling 
equations  may  be  reduced  to  one  equation.  This  is  done  by 
multiplying  the  mass  balance  by  B  and  subtracting  it  from 
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(9) 


(10) 


(11a) 


(lib) 


•  3k<  ■ 


f- 


the  energy  balance  to  give  an  equation  whose  only  solution  is 
X  2  =  B  x  x  (12) 

This  relationship  makes  the  reaction  rate  a  function  r>*  only 
one  variable. 

The  parameters  each  have  specific  physical  meaning.  The 
quantity  B  is  a  dimensionless  adiabatic  temperature  rise  and 
Da  represents  the  ratio  of  reactor  space  time  to  the 
characteristic  reaction  time.  Pe^  and  Pe^  are  the  Peclet 
numbers  for  mass  and  heat  transport  and  B  is  a  dimensionless 
heat  transfer  coefficient.  As  mentioned,  Le  is  the  ratio 
of  the  physical  transport  thermal  time  constant  to  the  physical 
transport  material  time  constant  167). 

Based  on  common  exothermic  reactions,  Hlavacek  and  Votruba 
[5,  Table  6.6]  list  values  of  B  and  Y  in  the  range  5-30  where 
the  high  values  of  B  usually  correspond  to  oxidation  reactions. 
For  the  "empty"  tubular  reactor  Le  =  1 ,  while  for  a  typical 
fixed  bed  reactor  Le  -  500.  When  the  reactor  is  empty, 

Pe^  =  p^2’  while  for  the  packed  reactor,  Pe^  ~  2-3  times 
Pe^  since  the  value  of  Pe^  derives  from  the  width  of  the 
reaction  zone  as  well  as  the  length  of  the  reactor.  This 
implies  that  multiplicity  is  possible  even  in  long  reactors 
(cf.  [41]).  The  Damkohler  number.  Da,  and  the  dimensionless 
heat  transfer  coefficient,  S,  vary  with  the  reactor  space 
time  but  are  usually  less  than  0.5  and  5  respectively. 

3.  BIFURCATION  ANALYSIS 

Bifurcation,  or  branching  of  solutions,  is  closely  related 
to  the  stability  and  thus  to  the  eigenvalues  of  linearized  system 


equations.  As  an  example,  consider  the  system  of  nonlinear 


coupled  ODEs 

dx 

di  = 

where  x  is  a  state  vector  and  p  represents  a  vector  of 
parameters.  The  system  is  locally  stable  if  all  eigenvalues 
of  the  Jacobian  have  negative  real  parts.  The  eigenvalues  are 
functions  of  the  system  parameters  and  these  may  change  such 
that  the  system  loses  its  stability.  The  exchange  of  stability 
occurs  as  some  eigenvalues  cross  the  imaginary  axis,  and  it  is 
at  this  point  that  the  bifurcation  can  take  place.  Therefore, 
bifurcation  is  often  referred  to  as  "the  principle  of  exchange 
of  (linearized)  stability"  (cf.  [68)  for  a  detailed  discussion 
of  this).  There  are  two  ways  in  which  the  eigenvalues  can  cross 
the  imaginary  axis,  namely: 

(i)  A  simple  eigenvalue  passes  through  the  origin.  This 
leads  to  bifurcation  of  stationary  solutions  and  is 
known  as  "static  bifurcation"  [69]. 

(ii)  A  pair  of  simple  complex  conjugate  eigenvalues  cross 
the  imaginary  axis.  This  leads  to  bifurcation  of 
periodic  solutions  and  is  known  as  " Hop f - b i f u r ca t i on " 
(cf.  (70)  for  more  details)  . 

Occasionally  multiple  eigenvalues  cross  the  imaginary  axis 
together  and  this  may  produce  complex  interactions  between  the 
two  basic  bifurcation  phenomena  (cf.  (70,71)  for  examples). 

Because  partial  differential  equations  (FDEs)  can  be 
regarded  in  some  aspects  as  an  infinite  set  of  ODEs,  one 
intuitively  expects,  and  can  in  fact  show,  under  certain 
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limited  conditions  (71-75  and  references  within] ,  that  the 


bifurcation  theorems  for  ODES  can  be  extended  to  the  "infinite 
dimensional  case".  The  review  paper  by  Crandall  [75]  gives  a 
particularly  readable  account  of  the  necessary  concepts  and 
proofs.  The  conditions  are  all  satisfied  for  the  parabolic  or 
elliptic  partial  differential  equations  which  arise  in  reaction- 
diffusion  problems. 

We  now  linearize  the  equations  around  a  spatially  varying 
steady  state  profile  rather  than  a  point  as  in  the  ODE  case. 

The  linearized  equations  take  the  form: 


3x 

3f(x) 

S  aF- 

Lx  +  ' 

--  dx 

x(z) 

(13) 

x(z)-x  (z) 

-  •  s 

B0xC°) 

«*  0 

(14) 

B.x(l) 

-  o 

(IS) 

-1- 

where  x  =  x-x  and  C,  L , 

-  -  s  - 

f,  Bq  ,  and  B  ^  are  defined  by 

equations  (8),  (9),  (10),  and  (11)  respectively.  The  branching 

to  nontrivial  solutions  or  periodic  solutions  is  then  governed 
by  the  discrete  spectrum,  i.e.  the  eigenvalues,  of  the  linear 
ope  rator  : 


F 

-x 


C 


-1 


L  + 


3f(x) 


3  x 


(16) 


Analogously  to  the  ODE  case,  one  sees: 


(i)  Static  bifurcation,  when  a  simple  eigenvalue  in 


the  spectrum  (’asses  through  the  origin. 

(li)  Hop f -bi f ur cat  ion ,  when  a  simple  pair  of  complex 
eigenvalues  in  the  spectrum  passes  the  imaginary 
axis. 

In  addition,  there  must  be  an  exchange  of  stability  at  the 
bifurcation  point;  e.q.,  in  order  to  have  Hopf-bi furcation , 
the  remaining  eigenvalues  must  have  negative  real  parts. 

The  first  N  eigenvalues  of  F  may  be  determined  by  using 
projection  techniques  such  as  Galerkin's  method  orthogonal 
collocation.  In  the  limit  of  infinite  terms,  both  methods 
will  represent  F^  exactly  (76,77],  Thus,  this  “late  lumping" 
procedure  retains  all  the  information  in  the  original  partial 
differential  equations,  contrary  to  the  behavior  of  the  “early 
lumping"  procedures  such  as  Hlavacek' s  “linearization"  [  3  B  j 
which,  for  example,  cannot  predict  the  existence  of  more  than 
three  steady  states.  For  both  the  collocation  ar.d  the 
Galerkin  method  the  finite  dimensional  Jacobian  takes  the  form; 
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where  N  is  the  number  of  eigenfunctions  in  the  Galerkin 
method  or  the  number  of  interior  collocation  points  in  the 
collocation  method. 

The  quantities  M.  and  K.  are  2x2  matrices  defined 
below  (except  in  the  special  adiabatic  case  with  Le  =  1, 
where  they  are  scalars).  In  the  following  we  shall  describe 
our  approach  to  this  analysis  using  weighted  residual  methods 
( i )  Galerkin 1 s  method 

We  make  the  usual  transformation  which  causes  the  problem 
to  become  self  anoint  [62]  : 

(  Pek  1 

Xk  ”  XkeXpi  '  ~2  Z  \  '  k  “  1  '  2 


and  choose  the  trial  function  expur. sion: 
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where  the  >  (*)  are  the  orthonor,,’a^  eigenfunctions 

corresponding  to  the  self  adjoint  eigenvalue  problem: 
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(77)  . 
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(ii)  The  Orthoqonal  Collocation  Method 


This  method  has  been  successfully  applied  to  many  chemical 
reaction  engineering  problems  similar  to  the  present  eigenvalue 
problem  (53,54*78-811.  In  this  method  the  first  and  second 
spatial  derivatives  are  approximated  by  a  weighted  sum  of  the 
values  of  the  dependent  variable  at  the  collocation  points: 


3x, 


3  z 


N+- 1 

l  A, 


j  =  0 


ij  kj 


i  N  +  l  ; 


1,2  (27) 


and 


32  x. 
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1,2  (28) 


The  weight  matrices  depend  on  the  trial  functions  which  in  our 

(a,  8) 


case  are  the  first  N  Jacobi  polynomials 


(z)  with  weight 


function  za(l-z)®  [81, ct.  3).  Since  there  is  no  special 


symmetry  in  the  tubular  reactor  problem,  we  use  (a, 8)  =  (0,0) 
in  which  case  the  polynomials  are  Legendre  polynomials.  In 
fact.  Georgakis  e_t  a_l-  154  ]  compared  the  convergence  rate  of 

the  eigenvalue  calculations  for  various  choices  of  a  and  8 
and  found  that  the  fewest  number  of  collocation  points  were 


re  j  ;ired  for  a  =  8  =  0.  The  collocation  points  are  then  the 

(0,0) 


zeros  of  P  '”'”'( z )  .  By  discretizing  the  equations  and 
N 


eliminating  the  boundary  conditions  (as  detailed  in  (54  ,81 

ct.  41),  we  obtain  the  following  set  of  ordinary  differential 

equations : 


where  the  eigenvalues,  X^n<  are  tbe  zeros  °f  the 
transcendental  equation: 
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By  inserting  the  trial  functions  into  the  linearized  equations 
and  making  the  resulting  residuals  orthogonal  to  the  first  N 
eigenvalues,  one  obtains  the  following  system  of  ODEs : 
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Note  that  j  =  5ji‘ 
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(ii)  The  Orthogonal  Collocation Method 


This  method  has  been  successfully  applied  to  many  chemical 
reaction  engineering  problems  similar  to  the  present  eigenvalue 
problem  153,54,78-81).  In  this  method  the  first  and  second 
spatial  derivatives  are  approximated  by  a  weighted  sum  of  the 
values  of  the  dependent  variable  at  the  collocation  points: 


and 


N  +  l 


j-0 


ij  kj 
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The  weight  matrices  depend  on  the  trial  functions  which  in  our 

case  are  the  first  N  Jacobi  polynomials  (z)  with  weight 

a  6 

function  z  (1-z)  [81, ct.  3).  Since  there  is  no  special 

symmetry  in  the  tubular  reactor  problem,  we  use  (a, 6)  =  (0,0) 
in  which  case  the  polynomials  are  Legendre  polynomials.  In 
fact,  Georgakis  e_t  a_l.  (54  1  compared  the  convergence  rate  of 
the  eigenvalue  calculations  for  various  choices  of  a  and  8 
and  found  that  the  fewest  number  of  collocation  points  were 
required  for  a  =  8  =  0.  The  collocation  points  are  then  the 
zeros  of  discretizing  the  equations  and 

eliminating  the  boundary  conditions  (as  detailed  in  [54  ,81 

ct.  4)),  we  obtain  the  following  set  of  ordinary  differential 


equations : 


where  y 


is  the  conversion 
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5. .  is  Kronecker's  delta  and  A  and  B  are  the  differentiation 

i}  - 

weight  matrices  corrected  for  the  boundary  conditions,  (14)  and 
(15),  i.e., 

^lj  "  Xij 

■I(A00'Pek)  Vl,N+l“A0fN+lVl,0]  1[  (AN+l.tH-lAOj"A0,N+lAN+l,J)Xi0 
+  ((Aoo-pek)Vi.j-Vi,oAoj)xi.JH.iI  (32) 


The  evaluation  of  the  derivatives. 


usually  requires 


3  f  . 


interpolation  among  the  discrete  points  of  the  steady  state 
solution.  The  computational  effort  involved  in  the  two 
procedures,  (i)  and  (ii),  are  then  nearly  equivalent  in  spite 
of  the  integrals  in  (26)  because  the  integrals  can  readily  be 
evaluated  by  quadrature.  The  choice  of  method  therefore 
strictly  depends  on  their  convergence  properties.  McGowin  and 
Perlmutter  [53]  showed  in  numerical  examples  that  Galerkin's 
method  converged  mono ton i ca 1 1 y  while  the  collocation  procedure 
converged  in  a  dampened  oscillatory  manner,  but  they  did  not 
compare  the  overall  rate  of  convergence.  Because  of  the 
similarity  between  the  e igen f unc tion  expansion  and  the 
perturbation  solution  to  the  special  adiabatic  case  with 
Le  =  1  [64,65]  ,  one  expects  that  the  necessary  number  of 
terms,  N ,  will  increase  with  the  value  of  the  Peclet  number 
starting  from  N  =  1  at  small  Peclet  numbers.  This  is  further 
discussed  below. 

4.  CALCULATION  OF  THE  STEADY  STATE 

The  steady  state  equations  for  the  reactor  (cf.  Equation  (7)) 
form  a  nonlinear  two  point  boundary  value  problem  which  is  quite 
stiff  even  for  moderate  values  of  Peclet  numbers  (Pe  ~  10)  . 

This  boundary  value  problem  must  be  solved  quickly  and  accurately 
if  one  wishes  to  obtain  the  bifurcation  curves  with  judicious  use 
of  computing  time.  Therefore,  we  have  given  special  consideration 
to  the  calculation  of  the  steady  state  profiles. 

Because  of  the  stiffness,  finite  differences  would  require 
far  too  many  mesh  points  and  can  thus  be  ruled  out.  Although 
efficient  routines  exist  for  solving  stiff  initial  value 
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problems  ( e .  g .  [82,83]),  t  h  >  ■  sheet.  i:«‘!  methods  of  McGinnis  (84) 

and  Kubicek  and  H  1  a  v  a  0  e  k  [  v  f>  ]  are  c  u  ;:iiu:  r  some  sir.  ce  they 
respectively  requite  the  integration  of  12  a  n d  24  first  order 
equations.  In  addition,  because  of  the  marching  nature  of  the 
technique,  the  solution  must,  be  stored  in  arrays  in  order  to 
be  used  in  the  bifurcation  analysis.  On  the  other  hand,  ir.  the 
method  of  weighted  residuals  the  solution  is  characterized 
by  a  few  trial  functions.  Moreover,  as  mentioned,  the  colloca¬ 
tion  approach  has  been  proven  to  be  very  efficient  in  solving 
this  type  of  boundary  value  problem.  There  are  several  ways 
this  method  may  be  applied  and  in  the  following  paragraphs 
we  review  the  advantages  and  disadvantages  of  these.  More 
details  on  the  computational  procedures  may  be  found  elsewhere  [37] 
(i)  Collocation  over  the  whole  domain,  0  5  z  £  1,  with  orthogonal 

polynomials.  This  method  is  simple  and  has  the  advantage  that 
the  collocation  points  and  weights  need  only  be  calculated  once. 

It  has  been  used  successfully  in  cases  of  moderate  values  of  the 
Peclet  numbers  (Pe  <  10)  [e.g.  54,79,80],  but  the  procedure  is 

inadequate  at  high  values  of  Damkohler  and  Peclet  numbers  where 
the  reaction  is  complete  within  a  narrow  zone  close  to  the 
reactor  inlet  as  illustrated  in  Figure  1.  In  order  to  obtain 
a  good  representation  of  the  narrow  reaction  zone,  a  large 
number  of  collocation  points  are  needed,  most  of  which  are  wasted 
downstream  from  the  reaction  zone.  Moreover,  the  large  number 
of  points  cause  the  interpolating  polynomial  to  wiggle  as  shown 
in  Figure  1.  Although  the  wiggles  are  slight,  they  significantly 
alter  the  bifurcation  behavior.  This  is  illustrated  in  Figure  2a, 
which  shows  the  calculated  regions  of  multiplicity.  Note  that 
the  collocation  technique  with  N  >  8  can  accurately  represent 
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AXIAL  POSITION 


Figure  1  Behavior  of  the  interpolating  polynomial  (a)  compared 
to  the  exact  solution  (b).  Approximate  solution 
based  on  orthogonal  collocation,  M  =  10.  Reactor 
parameters:  B  =  10.0,  Da  =  0.13,  Pe  -  Pe  =  15.0, 

B=  0.0,  Y  =  20.0. 


the  branch  corresponding  to  the  lower  steady  state,  but  even 


for  :i  -  10,  it  cannot  predict  the  steep  upper  steady  state 
profiles  when  Pe  >  10. 

( i i )  Transformation  of  the  independent  variable  followed  by 
collocation  over  the  whole  domain.  In  order  to  avoid  using  a 
large  number  of  collocation  points  when  the  profiles  are  steep, 
a  suitable  transformation  of  the  independent  variable  can  be 
made  such  that  the  reaction  zone  is  stretched  and  the  "dead- 
zone"  is  compressed.  This  implies  that  the  majority  of  the 
collocation  points  will  be  lo-ated  in  the  region  of  rapid 
changes.  However,  no  general  transformation  will  fit  the 
entire  region  transversed  by  the  reaction  front  for  changes 

in  the  bifurcation  parameters,  e.g.  Da.  The  transformation 
must  contain  at  least  one  adjustable  parameter  which  has  to 
be  fitted  to  an  approximate  profile  by  nonlinear  least  squares. 
Our  experience  with  this  method  indicates  that  the  computations 
become  too  lengthy  to  act  as  a  basis  for  the  bifurcation 
analysis. 

(iii)  Orthogonal  collocation  with  exponentials.  Orthogonal 

a.  8  .  z 

collocation  with  trial  functions  of  the  form:  f^(z)  -  z  e 

have  been  shown  to  give  excellent  results  in  plug  flow  reactor 
problems  since  the  trial  functions  are  solutions  to  the  corres¬ 
ponding  linear  problem  [86].  Because  the  axial  dispersion 
model  in  the  limit  of  large  Peclet  numbers  approaches  the  plug 
flow  model,  this  approach  should  also  be  able  to  solve  the 
stiff  cases.  However,  we  did  not  find  any  significant  improve¬ 
ment  over  the  standard  procedure  (i),  presumably  because  of 
difficulties  in  determining  an  optimal  choice  for  the 
coefficients,  and  8^.  The  optima' 


lection  of  these 


coefficients  remains  an  open  problem  [06) .  The  method  has  the 


further  disadvantage  that  re-evaluation  of  the  collocation  points 
and  the  weight  matrices  is  necessary  after  eacii  change  in  the 
bifurcation  parameter,  which  greatly  increases  computation 
time. 

( i v )  Orthogonal  Spline  Collocation  (also  referred  to  as 
orthogonal  collocation  of  finite  elements) .  This  method  has 
been  successfully  applied  in  other  problems  with  steep 
gradients,  notably  the  catalyst  particle  problem  for  large 
values  of  the  Thiele  modulus  [81,  ct.  7 ,  87,88).  The  increased 

accuracy  of  this  approach  over  the  global  approach  ti)  derives 
from  the  concentration  of  the  collocation  points  in  the 
regions  where  the  gradients  are  steep.  Since  the  number  of 
points  are  reduced  in  the  remaining  regions  where  the  solutions 
only  change  slightly,  the  total  number  of  collocation  points 
may  be  less  than  required  in  the  standard  case  (i).  In 
addition,  by  balancing  the  number  of  collocation  points  and 
spline  intervals,  one  can  obtain  a  patched  interpolating 
polynomial  of  sufficiently  low  order  that  wiggles  are  avoided. 

The  location  of  the  spline  points  is  clearly  critical  to  the 
accuracy  of  the  approach.  Carey  and  Finlayson  [87]  suggested 
placing  the  elements  such  that  the  mean  squared  residual  was 
minimized.  After  a  given  calculation  the  residuals  were 
examined  and  new  elements  inserted  where  the  residual  had 
been  largest.  Then  the  procedure  was  repeated  until  the  desired 
accuracy  was  reached.  This  approach  is  general  but  may  give 
large  array  sizes.  Instead,  one  may  fix  the  number  of  spline 
points  and  let  their  location  move  with  the  reaction  front  as 
the  bifurcation  parameter  changes.  A  simple  way  to  accomplish 
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this  is  to  monitor  the  gradients  of  the  solution  and  place  the 


spline  points  accordingly.  However,  this  approach  requires  a 
priori  knowledge  of  the  shape  of  the  profiles.  A  far  more 
general  approach  is  to  rearrange  the  spline  points  such  that 
the  mean  squared  residual  is  minimized.  This  can  be  done  by 
determining  the  sensitivity  of  the  solution  to  the  location  of 
the  spline  point  [37!.  Because  of  the  extra  computations  in 
the  op t i mi za t ion  of  the  spline  point,  the  choice  between  this 
method  and  the  one  proposed  by  Carey  and  Finlayson  depends  on 
the  relative  importance  of  computer  time  and  space  requirements. 
However,  for  cases  where  the  rate  expression  can  be  expressed 
as  a  function  of  a  single  reactant,  the  following  simple 
procedure  is  attractive. 

(v)  A  Simple  Orthogonal  Spline  Collocation  Method.  This  method 
only  applies  to  problems  where  the  reaction  rate  can  be  expressed 
in  terns  of  one  reactant,  but  because  of  its  ability  to  give  fast 
and  accurate  solutions,  it  is  worthwhile  considering.  Moreover, 
many  reactor  studies  involve  such  kinetics.  The  numerical 
difficulties  are  circumvented  here  by  obtaining  a  collocation 
solution  for  the  reaction  zone  and  patching  this  together  with 
an  analytical  solution  for  the  remaining  part  of  the  reactor, 
the  "dead-zone".  Thus  this  procedure  retains  the  advantages  of 
the  simple  collocation  procedures  described  in  (i)  above  and  in 
addition  provides  for  concentrating  the  collocation  points  in 
the  reaction  zone. 

If  one  assumes  that  the  reaction  is  essentially  complete 
beyond  the  spline  point  z  ^  ,  one  may  neglect  the  reaction 
rate  term  so  that  the  steady  state  equations  on  the  interval 
z  izfl  become: 

-i€>- 


with  the  boundary  condition 


z  =  1 


0 


and  the  continuity  condition 


(35) 


X 


2 


2s 


(36) 


where  x^s  is  specified  once  the  solutions  over  the  two  zones 
are  patched  together.  For  constant  x2w'  equation  (34b)  may  be 
solved  analytically: 


x  -x  „ 

2  2w 


"I6 


-a 2  ( 1- z ) 


-  a2e 


-o1(l-z) 


X2s'X2w  “lC 


-a  2 ( 1 -zs ) 


°2° 


-0l(l-zs) 
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whe  re 


Pe. 


a  .  a 
1*2 


1  ±  / l  + 
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Pe. 


(38) 


With  the  variable  change  Z,  =  z/z s>  the  steady  state 


equations  governing  the  reaction  zone  take  the  form: 


(39) 


with  boundary  conditions  : 


d5 


c-o 


*spVk 


k-1.2 


and  the  continuity  restrictions: 


(40) 


(41a) 


(41b) 


0  for  B=0 


(41c) 
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Following  the  sLb.-dard  collocation  technique  and  eliminating 
the  continuity  and  boundary  conditions,  one  obtains  the  colloca¬ 
tion  equations: 


M 

t 

j-1 


+  zs[fk(xliX2i)-62kBx2i^° 


1=1,..., M,  k=l,2 


(42) 
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where 


Lkij  zsAiJ 

+  [  (A00"Pck2s)  (\h-1  ,>H-l-bk0) 

BiO 

{lA0,JH-l(A>H-l,j"akO)  _  (AIH-l,fH-rbkO)AOj](pi^)  "  ZsAiO) 

M+l 

+  ^H-l,OA0j+(A0O+?ekZs^akO'AH<-l,j^  (  Pe^  zsAi,Mhl)} 


(43a) 
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for 


g>0 


(43d) 


For  a  given  set  of  parameters,  B,  Da,  Pe^,  Pe2  '  an<*  ^ 

these  equations  can  be  solved  by  Newton- Raphson 1  s  method  and 

the  spline  point  can  be  adjusted  such  that  x,  =  1  and 

1  ,  M+  1 

wiggles  are  avoided. 

In  order  for  the  method  to  handle  cases  where  the  reaction 
front  is  sharp  and  at  the  same  time  situated  close  to  the 
outlet  of  the  reactor,  one  or  two  additional  collocation 
points  may  be  needed.  Such  profiles,  of  which  Figure  3  gives 
an  example,  arise  when  the  Peclet  numbers  and  the  heat  transfer 
are  large.  The  steady  state  collocation  equations  are  then: 
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at  the  left  hand  boundary: 

M^l 

af0  "  AZl*klj 

at  the  interior  points  of  the  2'th  element: 

Mt+1  B 

*Q  [Y^  '  AziAllj)xklj  +  AzlIfk<xlli*  x2ti>  “  52kSx2ii 

i-1 . M  ;  1-1,. 

at  each  division  between  the  spline  intervals: 

*k, 1+1,1  : 

M.+l 

-L.  r  A  v  .  _1_ 

A*t  jIq  A2t+1 

l-l . L-l 


Tci.tt+1 


Mi+1+1 

j!0  Al+l.ljxk,Mfl,J' 


at  the  right  hand  boundary: 


AXIAL  POSITION 

Figure  4  An  example  of  steep  conversion  and  temperature 

profiles  found  by  spline  collocation.  (a)  Upper 
steady  state,  (b)  lower  steady  state.  B  =  16.8, 

Da  =  0.330,  Pc  *  320,  Pe  =  100,  x 


9 . 6 


FEED  TEMPERATURE  'C 


Ignition  and  extinction  behavior  of  outlet  conversion 
and  temperature  for  Lubeck’s  example  [69]  Pe  y  *  3 

Pe  j  *  100,  6  =  0.72  ,  B  =  8.7-106T^J,  Da  =  7 .  3 3 • 1 0  * e ~ ? , 

Y  *>  1.22  10I,T'1,  x,  =  yOlO-TlT;1  ,  T  feed  temperature 


as  the  feed  temperature  is  chunked.  Luieck  simulated  ignition 
when  the  feed  temperature  was  increased  from  4  P ') Q  C  to  4  6  4  °  C 
and  found  extinction  when  the  temperature  was  decreased  from 
3S5°C  to  330°C.  Our  results  are  in  agreement  with  those  values 
we  find  ignition  and  extinction  when  the  feed  temperature  is 
4  6  1. 0°C  and  3  3d. 2  0  C  respectively. 

The  computations  illustrated  in  Figures  3,  4,  and  5  were 

based  on  three  splice  intervals  (besides  the  dead-zone)  with 
six  collocation  points  in  each.  The  two  additional  spline 
points  were  not  really  necessary  in  the  calculation  of  the 
profiles  shown  in  Figure  4,  but  they  were  required  for  the 
calculation  of  the  intermediate  and  unstable  profiles  similar 
to  the  one  shown  in  Figure  3. 

5.  STATIC  BIFURCATION 

Three  methods  are  available  in  calculating  the  static 
bifurcation  points.  These  are:  (i)  a  direct  approach  entirely 

based  on  the  steady  state  collocation  equations,  (ii)  a  more 
general  extrapolation  procedure  based  on  finding  zero  eigen¬ 
values  of  the  Jacobian,  Equation  (17),  and  (iii)  a  turning 
point  calculation  based  on  the  Newton-Euler  Steady  State 
algorithm. 

( i )  The  direct  approach .  Here  the  bifurcation  points  are 
located  by  finding  the  parameter  set  for  which  the  Jacobian  of 
the  steady  state  collocation  equations  becomes  singular.  The 
bifurcation  point  can  only  be  approached  from  one  side  (in  the 
case  sketched  in  Figure  6  from  the  right  hand  side  on  the  upper 
branch  and  from  the  left  hand  side  on  the  lower  branch),  other¬ 
wise  the  Newton- Raphson  iteration  fails  to  converge  or 


converges  to  a  different  solution.  Therefore,  this  method  is 
not  recommended. 

( i i )  A  general  extrapolation  approach .  The  idea  here  is  to 

linearize  around  a  steady  state  uoint ,  (Da  ,  x  (z),  (z))  away 

s  Is  2  s 

* 

from  the  bifurcation  point.  Da ^ ,  and  then  use  the  condition 

det  (J)  =  0  to  successfully  extrapolate  to  the  bifurcation  value 

of  the  parameter  Da  (cf.  Figure  6).  This  requires  additional 

computation  in  evaluating  the  necessary  terms,  N,  in  the 

Jacobian  Equation  (17).  Figure  7  shows  for  various  values  of 

the  Damkohler  number  used  in  the  linearization,  Da  ,  the  value 

s 

of  the  predicted  bifurcation  point  DaL  as  functions  of  the 

b 

number  of  terms,  N,  in  the  Galerkin  expansion.  Note  that  while 
there  is  movement  toward  the  correct  value  of  Da,  iterative 
relinearization  and  prediction  are  required  for  convergence. 

However,  this  is  plagued  with  singular  problems  because 

* 

Da,  >  Da,  (cf  Fig.  6).  Expanding  the  dimension  to  incl’de  arc- 

D  D 

length  as  Keller  [90]  suggests  would  be  helpful,  but  we  four, 
the  following  method  to  be  even  more  efficient. 

(iii)  Turning  point  calculations.  Because  the  static  bifurca¬ 
tion  points  in  this  problem  and  similar  chemical  reaction 
engineering  problems,  e.g.  the  catalyst  particle  problem,  are 
turning  points,  one  may  use  a  very  simple  approach  involving 
suitable  change  of  dependent  variables  such  that  the  system 
equations  have  a  unique  solution.  This  approach  was  taken  by 
Sorensen  e_t  a_l.  [91]  who  were  able  to  trace  out  the  solutions 
to  the  catalyst  particle  problem  in  the  region  of  multiple 
steady  states  by  using  the  value  of  the  concentration  at  a 
collocation  point  as  a  parameter  rather  than  the  Thiele  modulus 
itself.  Here  we  also  choose  to  fix  a  concentration  and  include 
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the  Damkohler  number  as  a  dependent  variable.  In  the  adiabatic 
case  it  can  be  shown  that  the  profiles  cannot  intersect  [92), 
so  the  parameter  change  renders  a  set  of  equations  which  have 
a  unique  solution.  However,  in  the  nonadiabatic  case,  both 
the  concentration  and  temperature  profiles  can  intersect  [45, 
Figure  11).  Extensive  calculations  [e.g.  40,44,45,93)  indicate, 
as  is  physically  expected,  that  for  constant  E,  Pe^ ,  Pe^ ,  0, 
and  y  at  the  exit  of  the  reactor 


where  Da  is  a  dependent  variable.  The  boundary  condition 

* 

at  z  =  0  is  now  included  so  ,  .  and  L, . .  take  the  form: 

M  +  l ,  g  li j 


"  ^Wl.J  “  (A00  "  Pelzs)  AHH,0A0j 


(53a) 
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(53b) 


'lij 


hi  _ 

Pei 


2sAij  "  (A00_Pelzs)  A0j 


lO 

Pc, 


Z  A.a 
s  10 


while  L2 i j  *s  defined  by  (43a). 

This  set  of  nonlinear  algebraic  equations  can  be  solved 
by  the  Newton- Raphson  procedure,  provided  a  sufficiently  good 
initial  guess  is  available.  This  initial  guess  can  be 
determined  as  follows.  We  differentiate  the  equations 
with  respect  to  x,  „ . ,  to  determine  the  sensitivity  of 

i  i  M+I 

x,  . ,  x_  .  ,  Da  to  a  step  in  outlet  conversion  x,  • 
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(54) 


(55) 


(56) 
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This  system  of  linear  equations  can  be  solved  for  a  step  j.  n 

outlet  conversion  x ,  ,  to  yield  the  next  initial  guess  for 

1  /  M  + 1 

Eqns.  (50-52). 

For  the  first  step,  ore  may  begin  from  the  known  trivial 
solution  to  the  steady  state  equations- 


Pa  «  0: 


x  *.  0 
lj 

*2j  -  0  if  8=0 


Tt  B  Y 

X2J  2w 


1-Pe 


-a2(1-Zj)  'al(1-Zj) 

a. e  -  a-e 

1 _  z _ 

2  -a  -a 

2  2  2  1 
-  a2e 


(57) 


j  -  0.1,. 


The  system  of  equations  (50-52)  and  (54-56)  form  a  Newton- 
Euler  algorithm,  where  they  respectively  function  as  corrector 
and  predictor.  Beceise  the  left  hand  of  the  Newton-Raphson 
correction  solution  to  Equations  (50-52)  is  the  same  as  the 
left  hand  side  of  Equations  (54-56)  ,  the  inverse  of  the 
Jacobian  needed  for  the  Euler  prediction  step  is  already 
available  (if  for  example  the  LU  decomposition  is  used). 

Thus  the  prediction  step,  equations  (54-56),  requires  only 
an  inexpensive  back  substitution.  The  above  formulae  have  been 
based  on  the  simple  collocation  method,  but  two  additional 
spline  points  can  readily  be  included  as  was  done  above  in 
Equations  (44-48  ).  The  Newton-Euler  algorithm  is  •••sed  through¬ 
out  the  following  bifurcation  studies. 
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In  order  to  check  this  static  bifurcation  algorithm  we 

considered  two  previously  published  calculations.  Using  the 

spline  collocation  procedure  described  above,  we  calculated 

the  bifurcation  curves  shown  in  Fig.  2b.  These  are  in  good 

„  v 

agreement  with  the  results  of  Hlavacek  and  Hofmann  [38b,  Fig. 

151  except  for  a  very  slight  difference  in  the  value  of  the 
Peclet  number  at  the  trifurcation  point  (i.e.  the  confluence 
of  the  bifurcation  curves).  We  obtain  Pe  =  31.88  whereas 
Hlavacek  and  Hofmann  report  a  value  of  Pe  =  31.05. 

The  second  example  is  taken  to  be  a  case  with  five  steady 
states  studied  earlier  [93,  Table  2],  Our  procedure  gives 
the  same  values  for  Da  at  bifurcation  as  the  GPM-technique 
of  Kubicek  and  Hlavacek  [93].  This  example  assumed  y 
It  is  interesting  to  note  that  for  the  same  parameters  and 
rather  large,  finite  y  (y  =  100),  further 

calculations  show  that  only  three  steady  states  are  possible 
(in  agreement  with  recent  results  of  Kapila  e_t  a)^.  [46]). 

6.  HOPF  BIFURCATION 

Two  basic  approaches  exist  for  evaluating  the  Hopf  bifurca¬ 
tion  points:  (i)  direct  methods  bas' d  on  actual  computations  of 

the  eigenvalues  or  (ii)  indirect  procedures  based  on  the  magnitudes 
of  the  coefficients  of  the  characteristic  polynomial.  In  the 
first  method  all  the  eigenvalues  are  calculated,  preferably  by  a 
QR  algorithm,  (cf.  (94)  for  the  original  description  and  [81, 
p.  421)  for  an  efficient  FORTRAN  routine)  and  the  most  positive 
real  part  of  the  complex  eigenvalues  is  chosen  as  a  residual. 

This  is  then  made  as  close  to  zero  as  is  wanted  by  varying  the 
bifurcation  parameter  [95].  Clearly,  the  search  is  abandoned 


if  the  real  part  of  any  of  the  other  eigenvalues  is  positive, 
(ii)  In  the  indirect  methods  the  Hopf  bifurcation  points  are 
determined  by  varying  the  parameters  to  satisfy  necessary 
conditions  formulated  in  terms  of  the  coefficients  of  the 
characteristic  polynomial: 


P  t  X ) 


x2h-s1x2n-1+s2x2n-2- 


•+S2N-2X2-S2H-1X+S2N 


(58) 


Here  N  is  the  number  of  terms  in  the  collocation  or  Galerkin 
approximation  and  the  coefficients,  S^,  are  the  sums  of  the 
principal  minors  of  the  Jacobian  Equation  (17),  [96,  ct.  4). 

Necessary  conditions  that  the  characteristic  polynomial  has 
two  purely  imaginary  roots: 


*1.2  “ 


U)>  0 


are  : 


(-1) NuN+ (-1) N  1S2uN  x+. • *+ (-1) SjN2u+S2N  -  0 


(59) 


(60) 


and 


(-1)  1S1uN  1+(-l)N  2S3<jN  2+ - +  (-1)  S2N_3U+S2N-1  “  0  <61> 


In  order  to  have  bifurcation  we  must  further  have  an  exchange 

of  stability.  The  stability  test  is  conveniently  made  using 

the  Rou th- Hurwi tz  criterion  as  well  as  the  necessary  condition 

2  n-  i 

that  the  coefficients  (-1)  must  be  positive. 

For  N  ”  1  and  N  =  2  the  necessary  conditions  for  Hopf 
bifurcation  may  be  reduced  to: 

for  N= 1 

Sj  =  0  (62a) 

S2  >  0  (62b) 
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for  N  =  2 


S7  -  S2  +  S4  *  0 

■  1 

(63a) 

<  0 

(63b) 

<  0 

(63c) 

>  0 

(63d) 

Recently,  Kubicek  (97)  presented  two  variants  of  the  above- 
procedure  where  he  combined  the  Hopf  conditions  and  steady 
state  equations  and  solved  the  entire  system  by  Newt on - 
Raphson  iteration.  However,  such  a  procedure  is  likely  to 
cause  difficulties  in  the  present  tubular  reactor  problem 
because  of  the  large  search  space  combined  with  local 
convergence  of  Newton's  method.  Instead  we  use  the  Newton- 
Euler  steady  state  algorithm  and  search  over  Da  for  a  solu¬ 
tion  to  equations  (60)  and  (61)  in  the  parameter  region 

2  N-  i 

where  the  necessary  condition,  (-1)  >0,  is  satisfied. 

When  a  solution  has  been  found,  we  use  the  Routh-Hurwi t z 
criterion  to  check  the  sign  of  the  real  part  of  the  remaining 
eigenvalues.  This  approach  works  well  for  N  5  2,  whereas 
for  larger  N  the  direct  approach  (i)  gives  faster  and  more 
accurate  results. 

The  direction  and  stability  of  the  bifurcating  orbits  can 
be  determined,  in  principle, by  applying  the  standard  formulae  for 
ODEs  [70,98]  to  the  collocation  or  Galerkin  approximations 
(cf.  Section  3)  since  these  represent  F^  correctly  in  the 
limit  of  large  N.  Alternatively  one  may  use  the  new  formula 
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of  :i  and  I'o  a  I  <•  1  •  '  i  .  H  <.*  w  *  v  »  r  ,  t  h  ■  •  :  '.  .1 :  ■ 1  :  1  *  /  ! 

pet  i'll,  solul  iu:k>  1  *.  r .  >  ?  juibili1  i  in  iota:*  n  *'  r  •:  . 

At  too  t  mi  iliibt-  Hot-:  b  :  f  u  r  c  a  t  1  .  ■  n  a  !  ;*r  1  trims  w  >  r  v 

develuiitd,  no  ;dU'jl  jlioiis  of  If  op  f  b  j  f  or  c  a  t  1  on  points  *  r. 

tubular  reactor  equations  h ad  been  carried  oJt  to  our  knowledge. 

Therefore,  we  bad  to  be  satisfied  w  *  t  h  com;  arn.q  t  b*-  rararefer 

values  at  bifurcation  points  with  those  used  in  sir.  i! at  loro 

of  limit  cvcles.  for  example,  Hlavacek  and  Hofaar.n  [  3  V I 

found  a  stable  limit  cycle  for  the  parameters  B  -  11.0, 

Pe,  =  Pe„  =  1.0,  Le  =  1.0,  x„  =  0.0,  8  =  2.C,  Da  «  0.200  and 
12  2  w 

Y  -»  ®>.  This  compares  well  with  our  result  that  stable  limit 
cycles  should  exist  at  least  between  the  two  Hopf  bifurcation 
points.  Da  =  0.153  and  Da  =  0.238. 

In  addition,  the  forthcoming  paper  by  Heinemann  and  Poore 
[99]  offers  an  aposteriori  check  of  our  algorithms,  since  they 
have  calculated  the  Hopf  bifurcation  points  corresponding  to 
the  original  simulations  of  Varma  and  Amundson  [45b].  Applica¬ 
tion  of  our  algorithm  gave  the  same  results  as  obtained  by 
Heinemann  and  Poore  [99]  for  these  examples  [37). 

The  convergence  properties  of  the  ualerkin  and  collocation 
procedures  were  studied  in  general  by  varying  the  number  of  terms 

in  the  evaluation  of  the  critical  Lewis  number,  Le  .  This  number 

c 

represents  the  maximum  value  of  the  Lewis  number  that  allows  Hopf 

bifurcation  and  thus  gives  a  good  measure  of  the  convergence 

properties  of  the  above  algorithms.  The  actual  calculation  of 

Le  is  detailed  below.  As  also  found  by  McGowin  and  Perlmutter 
c 

[53],  we  show  that  the  Galerkin  procedure  converges  mono to n i c a  1 1 y 
while  the  collocation  method  converges  in  a  dampened  oscillatory 
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peripheral  sketches  in  Figure  J.  Up pal  et  a  1  .  [66]  also  give 

the  direction  arul  stability  oi  the  Li  furcating  orbits  as  well 
as  phase  portraits  of  the  possible  dynamic  behavior.  However, 
similar  results  for  the  tubular  reactor  are  not  included  here, 
because  classifying  the  stability  of  periodic  solutions  is 
beyond  the  sc  opt*  of  the  present  study. 

The  B-  8  parameter  ace  ;  lot,  illustrated  in  Figure  10, 
shows  the  var  us  regimes  of  bifurcation  behavior  for  the  axial 
dispersion  rr .  d  e 1  with  Pe  =  P e  ^  *  6.  This  diagram  was 
calculated  using  the  N e w t o n -  F u  1  o r  steady  state  algorithm  with 
2  spline  points  and  6-10  col  location  points  in  each  element. 

The  eigenvalue  calculations  were  based  on  Galerkin's  method 
with  the  first  6  eigenfunctions. 

In  addition  to  the  M 1 ,  S  ^  ,  and  S  ^  curves  also  found  in 
the  CSTR  case,  there  is  a  multiplicity  curve,  M ^  and  Hopf- 
bifurcation  curves  S^,  S  ^ ,  and  S  ^  .  The  curve  M  ^  represents 
the  appearance  of  a  second  set  of  static  bifurcation  points 
not  found  in  the  CSTR.  Note  that  for  large  3,  these  r.av 
occur  for  lower  values  of  B  than  the  first  set  (which  arise 
above  M^).  Thus  the  tubular  reactor  may  more  easily  show 
multiplicity  than  the  CSTR  for  large  amounts  of  cooling. 

The  curves  S^,  S  ^  allow  additional  Hopf  bifurcation 
points  to  appear  as  8  increases,  while  curve  S^  marks  the 
simultaneous  disappearance  of  two  Hopf  bifurcation  points  by 
coalescence.  These  M  and  S  curves  divide  the  parameter 
space  into  14  regions,  I-XIV;  the  first  6  of  which,  i.e., 

I-VI,  are  the  same  as  those  found  in  the  CSTR  analysis  (cf . 
Figure  9)  whereas  the  others  are  new.  The  bifurcation  behavior 
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i sners ion  mode 


expected  in  each  real  on  is  illustrated  by  the  sketches  of  the 
variation  of  the  steady  state  exit  conversion,  x^  ,  with  the 
bifurcation  parameter.  Da.  Specific  numerical  simulations 
of  changes  in  the  steady  state  exit  conversion,  x^,  and 
temperature,  x^,  with  Da  are  shown  in  Figure  12a-r  and  the 
corresponding  values  of  B  and  8  are  summarized  in  Table  2  and 
Figure  11.  The  vertical  lines  in  Figures  )2b,d,i  and  q 
indicate  the  values  of  Da  for  which  the  steady  state 
profiles  are  shown  in  Figure  13.  The  bifurcation  structure 
is  detailed  in  the  following  paragraphs. 

The  two  multiplicity  curves,  and  ,  intersect  at 

the  point  F,  (B,3)  =  (12.4,  2.2)  (cf .  Figure  10).  For  parameter 

combinations  to  the  left  of  P,  above  M ^  and  below  ,  x^ 

versus  Da  has  two  turning  points,  m2  ,  m^  ,  i.e.,  three 
steady  state  profiles  exist  for  a  range  of  Damkohler  numbers 
(cf.  Figure  12b).  Figure  1 3  ( i  )  gives  an  example  of  such 
multiple  steady  state  profiles.  If  (B,S)  moves  above  M ^  and 
,  a  second  set  of  static  bifurcation  points,  ,  m^  ,  appear 

in  the  ,  x^  curves.  In  this  case  five  steady  states  are 
possible  if  the  two  sets  of  multiplicities  overlap,  i.e.,  Da(m^) 
Da(rrij)  (cf.  Figure  1 2  i )  ;  otherwise  there  are  two  regions  of  three 

steady  states  as  illustrated  in  Figure  12p.  The  broken  curve 

* 

represents  the  parameter  combinations,  (B,S)  for  which 

Da(m^)  =  Da(m^)  and  thus  separates  the  regions  where  the  two 

* 

phenomena  occur.  M  is  not  a  bifurcation  boundary  but  it  is 
interesting  to  note  that  comes  very  close  to  the  point  0 

where  the  ,  S^,  anu  curves  intersect  (cf.  Figure 

To  the  right  of  the  intersection  of  and  M ^ ,  (point  P)  , 

-1.1- 


10)  . 


Table  2  Summary  of  numerical  simulations 
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three  steady  states  are  possible  for  ( B  ,  8 )  in  the  region  b e 1 o 
M  ^  and  above  M  ^ .  Thus  the  turning  points  are:  m  and  m  ^  , 
and  the  multiplicity  occurs  at  higher  conversions  than  is  the 
case  for  parameters  to  the  left  of  P.  Furthermore,  the 
multiplicity  now  occurs  to  the  right  of  the  maxi  mum  in  the  x^ 
curve  (compare  Figures  12b  and  d).  As  illustrated  by 
Figure  13(ii)  the  steady  state  temperature  profile  now  has 
a  maximum  inside  the  reactor  because  of  the  increased  heat 
transfer  (6).  This  maximum  moves  forward  and  becomes  more 
pronounced  as  the  reaction  ignites.  A  second  set  of  turning 
points,  and  m^ ,  appears  when  (B,B)  rises  above 

and  then  the  multiplicity  pattern  is  again  either  1-3-1-3-1 
(cf.  Figure  12q)  or  1-3  —  5-3  —  1  (cf.  Figure  12k).  In  the  very 
special  cases  where  Dafm^)  =  Da(m^)  or  Da (m^)  =  Da (m^)  the 
patterns  become  1-3-5-1  or  1-5-3-1,  respectively. 

Clearly,  if  (B,B)  pass  through  the  point,  F,  ,  m^,  m ^ , 
and  m^  appear  simultaneously  so  that  one  goes  from  a  region 
with  a  unique  steady  state  (cf.  Figure  12a)  to  a  region  where 
three  steady  states  exist  over  two  disjoint  intervals  of 
Damkohler  numbers  (cf.  Figure  12n).  Figures  13iii  and  13iv 
give  respectively,  examples  of  steady  state  profiles  for  cases 
with  1  —  3— 1-3  —  1  and  1  —  3  —  5— 3  —  1  multiplicity.  As  expected,  these 
profiles  combine  the  characteristics  of  the  ones  in  Figures 
13 (i )  and  1 3 (ii) . 

The  Hopf-bi furcation  structure  is  conveniently  denonstratet 
by  considering  the  different  phenomena  that  comes  about  as 
either  B  cr  8  is  varied.  First  1 ».  t  us  vary  B  keeping 
8  =  1.5.  In  region  I  we  have  a  unique  globally  stable  steady 
state  profile.  As  discussed  above,  two  turning  points  exist 


in  region  II  (cf .  Figure*  12b)  .  If  B  is  inert- a  sod  \  a:.t  if..- 
:>  1  curvv'  into  region  III  a  Hop  f  bifurcation  point,  s  ^  , 
slides  onto  t  h  o  upper  branch  (cf.  Figure  12c).  Above  , 

in  region  IX,  the  turning  points,  m  ^  ,  m  ^  ,  appear  to  the 
right  of  the  Hopf  point  and  therefore  do  not  interfere  with 
it  (cf.  Figure  1 2  j )  .  At  the  cur  '0  a  second  Hopf  point, 

i»  ,  ma  t  e  r  1  a  1  i  s  on  the  middle  branch  (region  XI,  cf.  Figure 
1  2 .:i )  ;  this  time  it  comes  in  from  m ^ . 

Since  region  XI  extends  to  the  upper  boundary  of  the 
parameter  space  plot ,  we  now  hold  B  constant  at  some  value 


:t  b  c  v  <  the  intersection  of  S  ^  and  S  ^  ,  say  15.5,  and  change  S. 
As  9  is  increased  the  Hopf  points  and  approach  each 

other  and  at  S4  they  coalesce  and  disappear.  This  means  that 
in  region  Vlllb  the  middle  branch  is  unstable  contrary  to  the 
'••ivc*  in  region  Villa  (compare  Figures  12h  and  1 2  i  )  .  Further 
2:1 -rouses  in  8  bring  us  past  where  a  new  Hopf  bifurcation 

pc.  1  nt  com*s  onto  the  upper  branch  from  (cf.  Figure  12k). 

Then,  8  passes  into  region  XIV  through  where  a  Hopf 

bifurcation  point  moves  past  m^  down  onto  the  lower  branch 
< f  Figure  12r).  Note  that  region  XIV  is  extremely  small  for 
*'•  4  .  Finally  increasing  B  leads  one  to  region  IV  through 
V. ,  ,  where  m^  and  m0  coalesce  and  disappear.  Figure  12d 
-  •  a  r*  m  p  1  i  f  i  e  3  the  structure. 

Next  let  us  fix  B  at  some  value  below  the  intersection  of 

?  an*.  5  but  above  M_  and  vary  B  again  starting  in 
>42 

re  r  i  r>  j*.  xi.  Now  S-  is  reached  before  SM  so  we  enter  region 
3  4 

/.ill  where  the  Hopf  point,  s^,  appears  on  the  upper  branch 
;  and  s  ^  still  exist  on  the  middle  branch  (cf. 


Figures  1 2  o  -  q )  .  A  further  increase  in  B  brings,  us  into 
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o  r 
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region  XIV  and  across  M  to  region  IV  as  before,  while  from 
roc  n  XII  we  enter  region  IV  directly.  If  M  ^  is  passed, 
the  multiplicity  vanishes  and  region  V  appears  (cf.  Figure  12e). 

Below  the  M_  curve  the  curve  forms  the  border  between 

regions  III  and  VI  where  a  second  Hopf  point  appears  on  the 

upper  branch  {cf.  Figure  12f).  If  B  is  increased  in  region  VI 
either  S  ^  or  M  ^  is  crossed.  In  the  first  case  (region  VII) 
a  third  Hopf  point  slides  onto  the  upper  branch  from  (cf. 

Figure  12g),  while  in  the  latter  case  (region  XII)  a 
multiplicity  is  formed  between  the  two  Hopf  points  (cf.  Figure 
12n).  Increasing  B  in  either  region  leads  to  region  XIII. 

Let  us  discuss  the  intersections  between  the  M  and  S 

•i 

curves.  Because  the  curves  and  S3  respectively 

represent  parameter  combinations  where  Hopf  bifurcation  points 
appear  from  under  the  turning  points  and  m^,  the  curves 

have  to  come  together  at  the  point  O  where  crosses  . 

Also  S ^  has  to  join  at  this  point  in  order  to  preserve  the 

I 

balance  of  Hopf  points.  Similarly  and  have  to 

coincide  where  S  ^  intersects  M  ^  since  they  represent 
points  which  come  onto  the  upper  branch  through  the  turning 
point  m.,.  Note  that  and  "kiss"  and  depart  again  as 


in  the  case  of  the  CSTR. 


A  rigorous  analysis  of  the  direction  and  stability  of  tit- 
branching  periodic  solutions  has  yet  to  be  made  in  order  10 
obtain  the  complete  phase  plane  structures  of  each  region. 

In  addition  to  the  nine  basic  phase  for  traits  given  by  l:  :  a] 
e  t  a  1  .  ([66],  Figure  6)  one  clearly  would  obtain  new  phase 

plots  in  the  regions  where  five  steady  profiles  and/or  three 
Hopf  bifurcations  exist.  However,  the  structure  of  these  w .  I  1 
not  be  entirely  new,  but  rather  combinations  of  characteristic 
of  the  basic  plots.  For  example,  the  phas»  portraits  of 
region  XIII  will  be  a  combination  of  those  for  regions  III  and 
IV.  Similarly,  the  plots  for  regions  IX  and  X  are  likely  to  b 
a  mixture  for  regions  II  and  III.  Because  a  unique  unstable 
state  is  surrounded  by  a  stable  limit  cycle,  one  knows  without 
further  analysis  that  stable  limit  cycles  exist  in  regions 
IV-VII  and  XII-XIV.  The  original  simulations  by  Varna  and 
Amundson  [45b,  Figures  2.3  and  2.4]  nicely  illustrate  with 
phase  portraits  and  reactor  profiles  some  of  the  possible 
limit  cycle  behavior  in  these  regions. 

Effect  of  Varying  Peclet  Number 

As  already  demonstrated  above,  the  most  dramatic  influence 
of  the  Peclet  number  on  multiplicity  and  oscillatory  phenomena 
is  the  introduction  of  an  additional  set  of  Hopf  and  static 
bifurcation  points  to  augment  those  found  for  the  CSTR. 
However,  it  is  also  rather  interesting  to  consider  the  effect 
of  the  Peclet  number  on  the  critical  values  of  B  and  B- 
Let  us  begin  with  the  multiplicity  limits  and  compare  these 
exact  limits  with  those  predicted  by  Hlavacek’s  "linearization 
[38]  and  the  one  point  collocation  method.  Both  these  early 
lumping  techniques  give  the  following  algebraic  steady  state 
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i?  xtfcnf.ivo  numerical  calculations  support  this  conjee  tut  o  . 

The  bifurcation  analysis  may  be  carried  through  with 
r  ct  to  Le  just  as  was  done  in  the  previous  sections  with 

respect  to  Da.  In  the  C  S  T  R  case  the  model  p a r  a m e  t e  r  s  can  be 
redefined  such  that  the  parameter  space  plots  for  Le  =  1  can 
be  used  directly  for  all  values  of  Le  167).  Unfortunately, 
in  the  tubular  reactor  case,  a  similar  redefinition  of  the 
parameters  cannot  be  made  because  of  the  convection  terms  i r. 
the  modelling  equations. 

When  N  -  1  in  the  collocation  or  Galerkin  approximation, 
it  is  possible  to  derive  an  expression  for  the  Lewis  number 
from  the  Hopf  condition,  Eqn.  (62a): 

Sx  *  tr(J)  =  tr(M^  +  Ku)  =  0  1691 

y  inserting  the  expressions  for  M  and  K  ^  (25,26)  and 

C  3  0 , 3 1 )  and  rearranging,  we  obtain  (for  the  Galerkin  procedure) 

the  expression 
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where  both  formulae  are  subject  to  the  constraint  Det(J)  >  0. 

When  N  -  2,  the  Hopf  condition  (63a)  may  be  rearranged 

to  a  fourth  order  Polynomial  in  Le  by  using  the  fact  that 

c 

the  Lewis  number  divides  all  elements  in  the  even  numbered 

rows  of  the  Jacobian  Equation  (17)  (371 .  However,  when  N  >  2, 

the  Hopf  bifurcation  condition  becomes  too  complex  to  be 

expressed  in  terms  of  Le  ,  and  Le  is  then  determined 

c  c 

iteratively  by  calculating  the  eigenvalues  of  the  Jacobian. 

In  the  following  paragraphs  we  give  examples  of  the  changes 
in  the  bifurcation  structure  as  Le  ■+  0.  Figures  16,  17,  and 

18  show  the  variation  in  the  critical  Lewis  number  with  the 
Damkohler  number  in  three  different  cases  of  steady  state 
behavior,  and  Table  3  summarizes  the  corresponding  bifurcation 
regions.  In  these  calculations  N  =  6  terms  were  used  to 
insure  accurate  values. 

In  the  first  casie  (Figure  16),  the  tubular  reactor  equations 

have  a  unique  solution  f or  all  values  of  Da.  This  is  globally 

* 

stable  when  L  e  >  L  e  where 
c 

Le 


max  L  c 
Da  >0 
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At  Le 


1.39  a  set  of  Hopf  bifurcation  points. 


3 


a  nd 


s  ,  appears  and  these  move  apart  as  I.e  decreases.  This 
behavior  corresponds  to  region  V.  Since  Le  has  a  second 

local  maximum,  Le  -  1.267  at  x,  -  0.72,  a  second  set  of 

c  1 

H o p  f  points  s  ,  and  s  e  m e  r  ci  e  s  w h  e  n  Le  -  1.267  .  Thus, 

12  c 

there  are  two  intervals  of  Da  numbers  w  1 1  e  r '  •  stable  o  s  c  i  1  1  a  t  i  o  r. 
exist.  This  structure,  which  was  not  found  above,  is 
classified  as  re  g  i_on__X_V  arid  is  illustrate  d  by  Figure  19b. 

When  Le  =  1.262  the  Hopf  points,  and  s^,  coalesce 

and  disappear.  Thus,  we  are  back  in  region  V. 

In  the  second  case  (Figure  17) ,  three  steady  state  profile 

* 

exist  for  0.160  <  Da  <  0.168.  As  Le  decreases  beneath  Le 

c 

Hopf  point,  appears  on  the  upper  branch  and  we  have  a  region 

III  behavior.  At  Le  =  1.27  a  Hopf  point  comes  onto  the 

lower  branch  at  the  turning  point  and  we  are  then  in  region  IV 

Becauso  Le  versus  Da  has  a  local  maximum,  Le  =  1.26  at 
c  c 

x.  =  0.71  and  a  local  minimum,  Le  =  1.22  at  x,  =  0.81, 

1  cl 

a  second  set  of  Hopf  points  emerge  on  the  lower  branch  when 

1.22  *  Le  5  1.26.  Thus,  there  are  now  three  bifurcation 

points  on  the  lower  branch  plus  one  on  the  upper  branch.  This 

structure  is  new  and  is  classified  as  region  XVI .  I t  i s 

illustrated  by  Figure  19c.  For  Le  <  1.22  the  set  of  points 

c 

disappear  again  and  we  are  back  in  region  IV. 

In  the  third  case  (Figure  18),  we  have  multiplicity  of  the 

* 

kind  1— 3-5-3-1.  Since  Le  is  located  on  the  middle  branch, 

c 

* 

wo  pass  from  region  VIII  to  XI  as  Le  decreases  beneath  I.e  . 

At  Le  =  1.27  the  Hopf  point  “falls  off"  the  second  turning 

point  and  then  we  have  region  IX  behavior.  At  Le  =  1.26  a  new 


b  i  f  u  r  c  a  1 1  o  n  point  omerqos  from  behind  the  third  turning  ;  >  j  : .  t 
arid  we  are  back,  in  t'vcjion  XI  .  When  be  =  1  .  2  2  a  Hwid  ;  ;  :it 

appears  on  tho  upper  branch  giving  us  region  XI  I  la  lic.idv  ior. 

At  be  -  1.20  the  Hopf  points  on  the  middle  of  the  hr  'inch 

c  o  a  I  e  s  c  e  and  we  h  a  ve  a  regio  n  X  a  s  tract  u  r  «,■  u  n  til  !.«•  -  'J  .  c  4  4 

where  a  point  emerges  on  tho  lower  branch,  (region  XIV)  .  "hi 
in  illustrated  by  Figure  1  9  a  ,  win  c  h  i  n  d  i  c  a  t  c  s  t  h.  •  *  •  x  i  s  f  •  r.  ■  •  ■  if 

a  stable  limit  cycle  surrounding  five  unstable  steal  y  s  t  j  t  < 
for  a  small  range  of  Da  just  to  the  riaht  of  the  lower  Ho;  f 
b i furcat ion  point. 

It  is  possible  to  answer  the  question  of  existence  of 
limit  cycles  in  fixed  bed  reactors#  by  considering  the  variati  • 
in  Le  with  Pe  as  illustrated  in  Figure  20.  Here-  also  the 

influence  of  the  number  of  terms  in  the  Galerkin  expansion  is 

* 

shown.  The  discontinuities  in  the  curves  occur  when  Le 

c 

shifts  among  branches  in  the  region  of  multiplicity  as  is 

* 

illustrated  by  Fiqures  16,  17  and  13.  Note  that  Le  becomes 

c 

less  than  unity  for  Pe  >  14#  which  implies  that,  even  for 
"empty"  tubular  reactors  (Le  =  1)#  oscillations  are  only 

possible  in  very  short  reactors  (as  was  shown  in  the  previous 
ctiord  .  Therefore#  oscillations  due  to  interaction  of  mixing 
and  reaction  terms  should  not  occur  in  industrial  fixed  bed 
:  -  rs  where  Le  “*  500  are  common . 

3  .  G 0 N C  r, l J D  1  N G_ R  EMA KK S 

I : .  *  : .  i s  p  a  p « •  r  wc  have  given  a  detailed  analysis  of  the 
•••sir.  ue  of  multiple  steady  states  and  oscillatory  behavior 
:  r.  tuinl.ir  reactors  using  the  pseudohomogeneous  axial 
Jisper  ion  model  a  •  an  example.  This  model  was  a  convenient 
choice  since  the  modelling  equations  formed  a  relatively  simple 


ft  of  non  1  inoar 


d  i  t  f  e  r  o  n  t  i  a  1 


coupled  p.trabol:  ■  partial 
fiju  it  !  orni.  Moreover  ,  p r  o  v  l  ou  s  1  y  p  uf  1  i  •_»  ho d  examples  of 
nu  It  1 1.'  1 '  s  t  o  a  d  y  starts  and  1  i  n  1 1  c  y  c  1  e  5;  f  u  »  t  h  i  s  mo  d  e  1  rr  i  J  e 
it  j'ossi  r  1  o  t  o  h  o  ck  t  ho  a  l  •  j  u  r  l  t  h  m  s  .  S;  »•  c  i  a  l  consider  a  t  i  =.  :• 

were  71  von  to  the  solution  of  stiff  steady  ■  tate  equations 
a  n  d  ,  a  no  n  g  d  i  f  f  e  r  o  n  t  approximation  roc  e  -.1  u  res,  orthogonal 
spline  collocation  was  found  to  give  the  1-st  accuracy. 
Comparisons  of  various  ways  of  calculating  bifurcation  to 
multiple  steady  states  showed  direct  turning  point  calculation 
to  be  most  convenient.  The  eigenvalues  of  the  linearized 
operator  were  determined  by  both  a  collocation  and  Galerkir. 
procedure.  The  former  converged  in  a  dampened  oscillatory 
manner,  while  the  latter  appeared  to  converge  faster  and 
mo  no t o  n i c  a  1 1 y  . 

The  bifurcation  algorithms  were  sensitive  to  the  accuracy 
of  the  steady  state  calculations  so  even  for  moderate  Peclet 
numbers  (Pe  **  5)  it  was  necessary  to  use  spline  intervals 
with  6-10  collocation  points  in  each  interval.  Although 
orthogonal  spline  collocation  was  very  efficient  in  the 
steady  state  calculations (  this  meant  that  10-60  minutes  of 
CPU  were  required  on  a  PDP  11/55  minicomputer  in  order  to  trac 
out  the  bifurcation  behavior  for  all  values  of  Da  >  0  w  h  c n 
all  other  parameters  were  fixed.  The  actual  length  of  t  h 
computations  depended  on  the  narrowness  of  the  search  interval 


a  nd 

the 

values  of 

B  and 

Pe  . 

To  fully 

do 

scribe  the  t  y  p »? 

d  yn  a 

n  i  c 

behavior  : 

possible , 

one 

still  has 

t  Q 

deternine  t  h e 

stab 

il.it 

y  and  dir* 

*ction  of 

the 

b  ranching 

O  f 

periodic  s 0 1 u  t 

This  can  easily  be  done  within  the  framework  of  the  present 


coefficient  in  Galerkin  expansion,  defined  in  Eq .  <10) 

collocation  weight,  see  Eq .  (27) 

dimensionless  adiabatic  temperature  rise,  defined  in 
Kq  •  (  0  ) 

c  -  •  1  location  we  i  g  h  t.  ,  see  Eq  .  (28) 

;  •'  n.idary  operator,  defined  in  Eq .  (11) 

ronccntrat ion 

heat  capacity  of  fluid  phase 
heat  capacity  of  solid  phase 
capacitance  matrix,  defined  in  Ea  .  (8) 

diameter  of  roactcr 

bamkohier  number,  defined  in  Eq .  (6) 

longitudinal  dispersion  coefficient 
activation  energy 

nonlinear  multivariable  function,  defined  in  Eq .  (10) 

linear  operator,  defined  ir.  Eq.  (16) 

-.r.thalpy  of  reaction 

lac obi  an, defined  Ln  Eq.  (17) 

pro -exponential  factor 

longitudinal  thermal  conduction 

?  x  2  natr ix»def ined  in  Eqs.  (26)  and  (31) 

reactor  length 

of  spline  points 

I  operator, defined  in  Eq.  (9) 

combined  collocation  weight,  see  Eq.  (43) 

I.ewis  number,  defined  in  Eq .  (6) 

• :  r  :  v  i  c.  a  1  L  ?  v;  is  number 


;  1 1 : 


a  tc 


i  ’  o  r  of  Jut  e  not  c  u  1  1  o  c  a  t  i  o  n  p  o  j  n  t.  s  i  n  steady 
s  o  1  u  t  i  o  n 

M.  2  x  2  natri  x,  tiffined  in  I-js.  ( 2  S )  and  (30) 

N  number  of  terms  in  collocation  or  Galerkin  approximation 


vector  of  parameter  values 


p  (  X ) 

characterist 

ic  polynomial,  Eq .  (58) 

Pel 

Fcclet  riuzu.e 

r  for  mass  dispersion. 

defined 

in  Eq  . 

(6) 

?e  2 

Peclct  numbe- 

r  for  heat  dispersion. 

defined 

in  Eq  . 

(6) 

R 

universal  gas  constant 

s  . 

1 

i'th  Hop  f  b l 

furcation  point. 

S  ^  sum  of  i ' t  h  principal  minor  of  Jacobian 

t  time 

t  dimensionless  time,  defined  in  Eq.  (6) 

T  tempera  ture 

wall  te  mp  e  r  a  t  u  r e 
Tq  feed  temperature 

Ut .  heat  transfer  coefficient 

v  ^  linear  gas  velocity 

x  state  vector 

x  ^  conversion, defined  in  Eq.  (6) 

x  ^  dimensionless  temperature,  defined  in  Eq.  (6) 

x^^  value  of  at  the  spline  point,  z^,  see  Eq.  (36) 

*2 w  dimensionless  wall  temperature 

y  state  vector 

z  longitudinal  coordinate 

z  dimensionless  longitudinal  coordinate,  defined  in  Eq.  (6) 

z  ^  i'th  spline  point 

7.  splint*  point  separating  the  reaction  zone  from  the 

"dead"  zone,  sec  Eq .  (36) 
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r. :i  example  of  steep  conversion  and  temperature 
profiles  found  by  spline  collocation.  (a)  t.'ppor 
steady  state,  (b)  lower  steady  state.  B  =  16.8, 
Da  0.310,  ?e.  =  320,  Pt-  =  100,  x =  <*.6, 

S  -  0.72,  Y  ~  16.9.  2  2W 


Ignition  and  extinction  behavior  of  outlet  conversion 
and  temperature  for  Lubeck's  example  •  Pfc  ■  “ 

Pe,  =  10  0,  6  =  0.72,  B  -  9.7-106?-  J,  Da  =  7  .  3  3  ■  1  '  ■■  ~ , 

Y  =  1-22  10-T-S  x2w  =  y  ( 31  0-T,  T"*  ,  T  f*«d  UEfir,,.r. 
in  °K. 

Steady  state  value  of  x  ^  (  z  )  for  varying  Da.  The  carv«. 
has  two  turning  points,  m^  and  . 


The  predicted  bifurcation  points,  Dau ,  for  various 
assumed  values.  Da  ,  and  varying  number  of  terms,  N, 
in  the  Galerkin  expansion. 

* 

The  critical  Lewis  number,  Le^,  f°r  varying  number  cf 
terms,  N ,  in  the  collocation  procedure  (x)  and  in 
Galerkin 's  method  ( o ) . 

Classification  of  the  dynamic  behavior  of  the  CSTR 
in  t  h  e  parameter  space  B  -  [66]. 
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1?.  KEY  fC’HDS  fConf  mue  on  reverse  side  if  necessarv  arid  identity  by  block  niwl-er  - 

Bifurcation,  part,  ial  differential  equations,  chemical  reactor, 
computat ional  methods 


2i  A  f,  S  T  P  A  2  T  (Continue  on  evfrjc  side  tl  necessnry  end  identify  ^y  block  number- 

Methods  for  studying  the  bifurcation  behavior  of  tubular  i eactei s 
been  developed.  This  involves  the  application  of  static  and  iiopf  bitun 
theory  for  PDF,’.-;  and  the  very  precise  determination  of  steady  state 
profiles.  Practical  computational  methods  for  carrying  out  this  analys 
iiscussed  in  some  detail.  For  tiie  special  case  of  a  !  list,  crier,  u  nv 
reaction  in  a  tubular  reactor  with  axial  dispersion,  the  bifurcation  h"1 
is  classified  and  summarized  in  parametei  space  plots.  In  partioulat  t! 
influence  of  the  lewis  and  Pec let  numbers  is  investigated.  It  is  shewn 
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.  r  :  !uc  ib  interaction  of  dispersion  and  reaction  ef  fects  should  not  B 

•  i  xr'!  !  <■  i  r  (Victors  and  moreover,  should  only  occur  in  very  short  I 

:  i.i !  a*  i  eactors.  The  parameter  study  not  only  brinys  together  *1 

•  :  l.  :.<\1  examples  of  multiple  and  periodic  solutions  but  also 
I'o-rt;'  undiscovered  wealth  of  bifurcation  structures.  Sixteen  of 
•  lire:,  which,  come  about  by  combinations  of  as  many  as  four 
...  t  o  multiple  steady  states  and  four  bifurcations  to  periodic 
,  i:e  illustrated  with  numerical  examples.  Although  the  analysis  is 

•  ;  ... ■udohomouenuous  axial  dispersion  model,  it  can  readily  be 

.<  ot net  reaction  diffusion  equations  such  as  the  general  two  phase 

f  ;  Vi-  i  be  1  reactors.  ;; 
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